import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用黑体显示中文
plt.rcParams['axes.unicode_minus'] = False    # 正常显示负号

# 物理常数
g = 9.8  # 重力加速度 (m/s^2)

# 系统参数
l = 1.0    # 摆长 (m)
m = 0.1    # 质量 (kg)

# 有阻尼单摆的微分方程
def pendulum_eq(y, t, k):
    """
    有阻尼单摆的微分方程
    y[0] = theta (角度)
    y[1] = dtheta/dt (角速度)
    k = 阻尼系数
    """
    theta, omega = y
    dtheta_dt = omega
    domega_dt = -(k/m) * omega - (g/l) * np.sin(theta)
    return [dtheta_dt, domega_dt]